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Abstract 

This  paper  proposes  a  novel  numerical  methodology,  termed  the  discrete  thermal  element  method  (DTEM),  for  the  effec¬ 
tive  modelling  of  heat  conduction  in  systems  comprising  a  large  number  of  circular  particles  in  2D  cases.  Based  on  an  exist¬ 
ing  analytical  integral  solution  for  the  temperature  distribution  over  a  circular  domain  subjected  to  the  Neumann 
boundary  condition,  a  linear  algebraic  system  of  thermal  conductivity  equations  for  each  particle  is  derived  in  terms  of 
the  average  temperatures  and  the  resultant  fluxes  at  the  contact  zones  with  its  neighboring  particles.  Thus,  each  particle 
is  treated  as  an  individual  element  with  the  number  of  (temperature)  unknowns  equal  to  the  number  of  particles  that  it 
is  in  contact  with.  The  element  thermal  conductivity  matrix  can  be  very  effectively  evaluated  and  is  entirely  dependent 
on  the  characteristics  of  the  contact  zones,  including  the  contact  positions  and  contact  angles.  This  new  element  shares 
the  same  form  and  properties  with  its  conventional  thermal  finite  element  counterpart.  In  particular,  the  whole  solution 
procedure  can  follow  exactly  the  same  steps  as  those  involved  in  the  finite  element  analysis.  Unlike  finite  elements  or  other 
modern  numerical  techniques,  however,  no  discretization  errors  are  involved  in  the  discrete  thermal  elements.  The  mod¬ 
elling  error  mainly  stems  from  the  assumption  made  about  the  heat  flux  distribution  within  the  contact  zones.  Based  on 
some  theoretical  work,  an  enhanced  version  is  suggested  to  improve  the  approximation.  The  numerical  assessment  against 
finite  element  results  indicates  that  the  basic  version  of  DTEM  can  achieve  a  very  reasonable  solution  accuracy,  while  the 
enhanced  version  further  improves  the  accuracy  to  a  high  level.  In  addition,  thermal  resistance  phenomena  between  the 
contact  zones  can  be  readily  incorporated  into  the  current  modelling  framework. 
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1.  Introduction 

Heat  transfer  in  particle  systems  can  be  encountered  in  many  engineering  and  scientific  applications,  for 
instance,  in  porous  media,  soil  and  geo-mechanics,  civil,  chemical  and  nuclear  engineering,  and  materials  pro¬ 
cessing.  Although  effective  numerical  thermal  analysis  techniques,  predominantly  finite  element  method  based 
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(FEM)  [1],  have  been  well  established  for  various  types  of  problems,  the  modelling  of  heat  transfer  in  systems 
comprising  a  large  number  of  particles  remains  a  computational  challenge,  particularly  when  coupled  with 
other  physical  phases  such  as  solid  and  fluid  media. 

In  general,  two  different  approaches  have  been  proposed  to  undertake  thermal  analysis  of  problems  involv¬ 
ing  particles.  When  a  particle  system  can  be  treated  as  a  porous  medium,  a  ‘smeared’  continuous  approach  is 
often  adopted  where  the  effective  macroscopic  properties  (assumed  homogenous)  of  the  system  is  generally 
determined  by  a  detailed  thermal  analysis  using  the  finite  element  method  (FEM),  for  instance,  by  taking  a 
small  scale  micro-structure  of  the  system  or  a  representative  volume  element.  The  most  attractive  feature  of 
such  an  approach  is  that  it  can  utilize  the  established  continuum  based  methods  to  achieve  a  high  computa¬ 
tional  efficiency  for  large  scale  problems.  Nevertheless  the  approach  cannot  easily  account  for  structural  vari¬ 
ations  at  both  spatial  and  particle  levels.  There  are  a  number  of  emerging  applications  where  the  consideration 
of  features  at  all  scales  may  be  essential.  Theoretically,  any  type  of  heat  transfer  problem  can  be  modelled  by 
FEM.  However,  FE  methods  would  be  extremely  inefficient  to  model  the  entire  particle  system  due  to  an 
excessively  large  scale  finite  element  model  resulting  from  the  discretizition  of  each  particle.  In  this  context, 
the  boundary  element  method  [2,3]  seems  to  be  a  better  option  since  only  the  boundary  of  a  particle  needs 
to  be  discretized,  which  gives  rise  to  a  much  smaller  numerical  system  in  comparison  with  the  FEM,  though 
this  possibility  has  not  been  fully  exploited,  to  our  best  knowledge.  Clearly  more  efficient  thermal  modelling 
technologies  are  required  for  large  scale  particle  systems. 

On  the  other  hand,  some  attempts  have  been  made  in  which  the  particles  are  individually  represented  but  by 
a  much  simplified  model.  One  notable  approach,  termed  the  thermal  network  [4],  represents  each  particle  as  an 
isothermal  disk/sphere  and  each  pair  of  contacting  particles  by  an  equivalent  thermal  resistance  bar  or  pipe  so 
that  the  thermal  equations  of  the  whole  system  can  be  modelled  in  terms  of  the  corresponding  thermal  resis¬ 
tance/conduction  network.  Although  simple  and  computationally  efficient,  the  isothermal  assumption  made  in 
the  model  over  simplifies  the  actual  thermal  behavior  of  particles  and  limits  its  applicability  for  transient  cases 
[4].  In  a  different  approach  proposed  in  [5],  the  particle  system  is  represented  as  a  Voronoi  polyhedron,  and  the 
effective  conductivity  of  the  packing  is  modelled  again  by  network-type  analyzes.  The  common  feature  of  the 
above  two  approaches  is  that  the  determination  of  the  thermal  properties/parameters  in  the  models  is  based 
more  on  an  ad  hoc  manner  than  on  a  rigorous  theoretical  foundation,  which  may  introduce  modelling  errors 
that  are  difficult  to  quantify  and  control. 

Furthermore,  heat  transfer  within  a  particulate  system  is  often  closely  coupled  with  particle  mechanical 
contact  and/or  fluid  flows.  Over  the  last  decade  or  so,  many  coupled  numerical  techniques  have  been  devel¬ 
oped  to  account  for  the  interactions  between  different  physical  phases.  Some  latest  developments  have  focused 
on  the  accurate  modelling  of  different  phases  and  their  coupling  in  particle  transport  problems  (see  for 
instance  [6]).  Coupling  of  discrete  elements  with  Navier-Stokes  equations  for  particle-fluid  interaction  prob¬ 
lems  has  been  extensively  used  for  the  modelling  of  fluidized  beds  [7,8].  A  simple  particle-particle  heat  transfer 
model,  analogous  to  the  standard  discrete  elements  for  solids,  has  also  been  introduced  [9]  to  account  for  ther¬ 
mal  effects.  This  particle  thermal  model  is  nevertheless  similar  to  the  thermal  network  approach  mentioned 
earlier  and  thus  suffers  from  poor  accuracy.  In  the  coupled  Lattice  Boltzmann  method  (LBM)  [10]  and  the 
discrete  element  method  (DEM)  [11]  for  particle-fluid  interaction  problems  [12-15],  (turbulent)  fluid  flows 
are  modelled  by  the  LBM  while  individual  solid  particles  are  described  by  the  DEM,  and  the  particle-fluid 
interaction  can  be  accurately  calculated  through  a  suitable  boundary  condition.  Although  extension  to  include 
heat  conduction  has  also  been  achieved,  currently  available  thermal  Lattice  Boltzmann  formulations  [16,17] 
are  not  efficient  for  simulating  heat  conduction  in  systems  comprising  a  large  number  of  particles. 

The  main  objective  of  this  paper  is  therefore  to  propose  a  novel  numerical  methodology  for  accurate  and 
effective  modelling  of  heat  conduction  in  individual  particles  as  well  as  particulate  systems  in  2D  cases.  Under 
the  assumptions  that  particles  are  circular  disks  and  conduction  is  the  dominant  heat  transfer  mechanism  in 
the  problem,  the  proposed  approach,  termed  the  discrete  thermal  element  method  (DTEM),  provides  a  very 
simple  but  accurate  description  of  heat  conduction  within  the  system.  Consequently,  the  approach  can  be 
employed  to  replace  or  improve  the  accuracy/efficiency  of  existing  particle  thermal  models  and  to  very  effec¬ 
tively  simulate  heat  conduction  in  systems  comprising  a  large  number  of  particles. 

The  development  of  the  DTEM  is  motivated  conceptually  by  the  boundary  element  method  (BEM) 
where  only  the  boundary  of  a  particle  needs  to  be  discretized.  Its  final  form  closely  resembles  finite  element 
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formulations.  Specifically,  by  utilizing  an  existing  analytical  integral  solution  for  the  temperature  distribution 
over  a  circular  domain  subjected  to  the  Neumann  boundary  condition  (not  included  in  [18]),  a  much  simpler 
linear  thermal  conductivity  equation  is  established  for  the  particle,  between  the  average  temperatures  and 
resultant  fluxes  at  the  contact  zones  of  a  particle  with  its  adjacent  particles.  Thus,  the  only  unknown  variables 
for  each  particle  are  the  average  temperatures  at  the  contact  zones,  which  typically  number  around  3  or  4. 

The  paper  is  organized  as  follows:  in  Section  2,  the  integral  solution  to  heat  conduction  in  a  circular  domain 
is  introduced.  Then  the  details  of  the  proposed  DTEM  are  described  in  Section  3.  Additional  issues  associated 
with  the  DTEM  are  discussed  in  Section  4.  Two  examples  are  provided  in  Section  5  to  assess  the  solution  accu¬ 
racy  of  the  method  against  FEM  and  to  demonstrate  its  effectiveness.  The  main  features  of  the  proposed  meth¬ 
odology  are  summarized  in  Section  6. 


2.  Analytical  solution  of  temperature  distribution  in  a  circular  domain 


Consider  a  2D  circular  domain  of  radius  R  with  an  isotropic  thermal  conductivity  k  and  subjected  to  a  pre¬ 
scribed  heat  flux  (Neumann)  boundary  condition  as  shown  in  Fig.  1.  Using  a  polar  coordinate  system  with  the 
origin  set  at  the  centre  of  the  domain  Q  =  {(r,  6)  :  0  ^  r  ^  R;  0^0^  2n},  the  temperature  distribution  in  the 
domain  is  governed  by  the  Laplace  equation  as 

J  kAT  =  0  in  Q 
(Kf  =  4(0)  °n  8*2 


(1) 


where  denotes  the  boundary  (circumference)  of  the  domain,  and  n  is  the  outward  normal  to  the  boundary. 
Note  that  the  flux  q(6)  inward  to  the  boundary  is  assumed  positive.  The  temperature  at  any  point  (r,  6)  <G  Q 
can  be  expressed  by  the  so-called  Dini’s  integral  [19]  as 

r2n 


T{r • e) = -  IL  L  « w  in  ['  -  cos(e  -  + © 

where  C  is  a  constant.  The  existence  of  this  solution  requires 

pin 

q{6)d8  =  0 


d  <f>  +  C 


j 


(2) 


(3) 


which  is  physically  equivalent  to  the  global  heat  flux  equilibrium  of  the  domain.  It  is  easy  to  verify  that  the 
temperature  at  the  centre  T0  =  T( 0, 0)  =  C.  Thus, 

r2n 


T{r' B) = - £-k l  in  1  - 2'r - « + © 


d (p  +  T0  (r,  6)  G  £2 


(4) 


Fig.  1.  Circular  domain  with  Neumann  boundary  condition. 
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Then  the  temperature  on  the  boundary  is  obtained  as 

Tc{6)  =  T(R,6)  =  -f-  f  q(4>)  ln[2  —  2cos(0  —  4>)}d(j)  +  T0  (0  <  6  <  2n) 

271K  J 0 

By  utilizing  the  condition  (3),  the  above  formula  can  be  rewritten  as 


T< (0)  =  -£  /%(*) 
7ZK  J o 


In 


.  e-(j) 

sin — - — 


d  <fi 


(5) 

(6) 


The  significance  of  this  solution  is  its  explicit  evaluation  nature.  Recall  that  only  a  Fredholm  integral  equation 
of  the  second  type  arises  in  the  BEM  [2,3],  which  then  has  to  be  solved  numerically.  Note  that  solution  (6)  is 
the  basis  for  the  development  of  the  discrete  thermal  element  approach  described  below. 

According  to  the  mean  value  theorem  for  Laplace’s  equation  [19],  the  maximum  and  minimum  values  of 
the  temperature  are  always  achieved  on  the  boundary  6£2.  In  addition,  the  temperature  at  the  centre  of  the 
domain,  T0 ,  is  the  average  value  along  any  circle  centered  at  the  origin  within  the  domain: 

T°  =  ^l  T{r,e)&e  (0  <r^R)  (7) 

Particularly, 

To  =  5n  C  (8) 

It  can  be  deduced  that  T0  is  also  the  average  temperature  over  the  whole  domain: 

To  =  ^  [  T(r,e)dQ  (9) 

FK  Jq 


3.  Discrete  thermal  element  method 

Consider  a  particle  assembly  subjected  to  a  thermal  boundary  condition,  as  shown  in  Fig.  2a,  in  which  each 
particle  may  be  in  contact  with  a  number  of  adjacent  particles,  and  heat  is  conducted  through  the  contact 
zones  between  the  particles.  Note  that  it  can  be  considered  that  the  particles  are  brought  into  contact  not  only 
under  mechanical  forces  but  also  being  'cemented’  or  bonded,  as  illustrated  in  Fig.  2b. 

Let  us  isolate  an  arbitrary  particle  that  is  in  contact  with  n  neighboring  particles.  (Unknown)  heat  flux  dis¬ 
tributions  are  imposed  over  the  n  isolated  contact  zones  on  the  boundary  of  the  particle,  and  the  remainder  of 


(a)  Assembly 


Fig.  2.  Heat  conduction  in  a  simple  particle  system. 
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the  particle  boundary  is  assumed  fully  insulated,  as  shown  in  Fig.  3.  Here  the  contact  zones  are  assumed  as 
circular  arcs,  although  they  may  be  flat  or  other  shapes  under  mechanical  contact  [20].  This  approximation 
may  be  acceptable  as  contact  normally  occurs  over  a  small  region.  In  addition,  each  contact  zone  is  described 
by  both  the  position,  in  terms  of  the  angle  6  of  its  middle  point,  and  the  (half)  contact  angle  a  that  determines 
the  contact  arc  length.  Generally,  the  position  angles  are  well  spaced  along  the  boundary  and  the  contact 
angles  a*  are  small.  For  a  dense  packing,  the  average  coordination  number  of  particles,  i.e.  the  value  of  n ,  is 
around  3^1.  The  positions  and  contact  angles,  0*  and  a h  of  all  the  n  contact  zones  constitute  the  local  element 
(thermal  contact)  configuration  of  the  particle.  Note  that  the  particle  radius  R  plays  no  direct  role  in  the  for¬ 
mulation,  as  will  be  shown  later. 

Next  the  heat  conduction  over  the  entire  particle  is  to  be  represented  with  a  discrete  relation  between  the 
average  temperatures  and  resultant  fluxes  at  the  n  contact  zones. 


3.1.  Discrete  temperature  distribution  on  particle  boundary 

The  boundary  flux  condition  in  (1)  can  be  refined  for  the  special  case  shown  in  Fig.  3  as  follows: 

’  qt(6  —  9i )  9i  —  a*  ^  0  ^  0*  +  a*  (/=!,...,  n) 


q(9)  = 


0 


otherwise 


From  (6),  the  temperature  on  the  boundary  becomes 

w  pO 


Tc(9)  =  —  — 

7LK 


n  pffi+OLi 

E  /  M) 

j=  1  JOj-aj 


In 


.  e-cj) 

sin — - — 


d  <f>  +  T0 


or 


/?  . n  n 

7=1  J  ~aJ 


In 


sin- 


d(/>  +  T0 


and  the  temperature  distribution  along  the  ith  contact  arc  is 

0  —  4>  —  0j 


n  n  pa 
TIK 


n  r*j 

E  / 

/= i j -y-i 


In 


sin- 


d(p  +  T 0  (0j  —  a,  ^  0  ^  0j  +  a,) 


(10) 


(11) 


(12) 


(13) 


Ot  +a; 

Q t=R  [  q(6)d6 

0-a, 


(resultant  flux) 


Fig.  3.  Problem  description  of  heat  conduction  in  a  circular  particle. 
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By  defining  the  resultant  flux  Qt  and  the  average  temperature  Tt  on  the  zth  arc,  respectively,  as 

Q,  =  R  r  q,(B) de 

J —a  i 


i  rdi+di 

-  /  Tim m 

J  Oi—OLi 


T  = _ 

!  2a, 

and  assuming  that  qt{9)  is  constant  over  the  zth  arc,  i.e.  qt{9)  =  qt,  the  resultant  0,  becomes 

Q,  =  2ctiRqi 

According  to  Eqs.  (13)  and  (16),  the  average  temperature  7)  over  the  arc  can  be  obtained  as 


Ti  2a, ■ 


Y  rOi+cti  r  n  potj 

:a>  Jo,-*,  kk  ^  J_ 


<7  ; l  n 


sin- 


d  <f> 


d  9  +  T0 


or 


^  =  E 

7=1 


Qt 


An  K(x} 


»JT 

i^j  J  -a,-  J -oij 


In 


sin- 


A  Gy  +  9  —  (j) 


d4>d9 


T 

1  o 


(14) 

(15) 

(16) 

(17) 

(18) 


where  A 9tj  =  0,  —  9j,  which  is  the  angle  ‘distance’  between  the  zth  and  /th  arcs.  The  above  equation  can  be  ex¬ 
pressed  compactly  as 


Ti  =  ^2  hijQj  +  T0  (z=  1,...,zj) 

]=  i 


where 


*„=-—! —  r  r 

AnKotiOCj  J  J 


ln 


.  A  0^  +  6  —  (j) 
sin — - — - - 


dc/)de 


(19) 


(20) 


Eq.  (19)  essentially  establishes  an  explicit ,  linear  algebraic  system  of  equations  for  the  average  temperature  on 
a  contact  arc  and  all  the  (resultant)  fluxes  acting  on  the  contact  zones,  which  provides  the  theoretical  founda¬ 
tion  for  the  development  of  the  discrete  thermal  element. 

The  double  integrals  in  formula  (20)  are  inter-changeable  and  the  integrand  is  always  negative,  therefore 

(21) 

(22) 


hij  =  hp  >  0 


Also  note  that  the  heat  flux  equilibrium  condition  (3)  changes  to 

n 


i=  1 


3.2.  Discrete  thermal  element  equations 

By  introducing  the  element  temperature  and  flux  vectors, 

Te  =  {Ti, . . . ,  r„}T  and  Qe  =  Q„Y 

and  the  element  thermal  resistance  matrix  He, 

He  =  N„x„  (23) 

as  well  as  a  n  x  1  vector  with  all  components  of  unity,  e  =  {1, . . . ,  1}T,  Eq.  (19)  can  be  expressed  in  matrix 
form  as 


Te  -  eT0  =  HeQ( 


(24) 
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This  is  the  element  heat  conduction  equation  in  terms  of  thermal  resistance.  The  equation  possesses  the  fol¬ 
lowing  features:  (1)  the  resistance  matrix  He  is  symmetric  (cf.  Eq.  (21));  (2)  the  temperatures  Te  are  decoupled 
and  thus  no  simultaneous  equations  need  to  be  solved  for  determining  the  temperature  when  the  fluxes  Qe  are 
known.  These  two  properties  are  a  distinct  advantage  over  the  conventional  boundary  element  formulation  in 
which  the  corresponding  matrix  H  is  generally  unsymmetric  and  the  nodal  temperatures  are  coupled. 
Though  not  mathematically  proved,  He  is  invertible.  Let 

Ke  =  He  1  (25) 

then  the  inverse  form  of  Eq.  (24)  reads 

Ke(Te-ero)  =  Qe  (26) 

Note  that  in  both  (24)  and  (26),  the  discrete  temperatures  Te  are  relative  values  to  the  average  temperature  T0 
which  is  as  yet  unknown.  Thus,  a  further  treatment  is  required  so  that  equations  (24)  and  (26)  are  fully 
solvable. 

The  total  heat  balance  condition  (3)  can  be  rewritten  as 

eTQe  =  0  (27) 

Multiplying  (26)  by  eT  gives 

T0  =  glTe/fce  (ge  =  Kee,  Ke  =  eTge)  (28) 

i.e.  the  average  temperature  of  the  particle  can  be  readily  obtained  by  a  linear  combination  of  the  discrete 
boundary  temperature  Te.  This  new  formula  can  thus  be  viewed  as  the  discrete  version  of  the  mean  value  con¬ 
dition  (28).  In  addition,  Ke  may  be  interpreted  as  the  total  element  conductivity. 

Based  on  relation  (28),  eliminating  T0  from  Eq.  (26)  leads  to  the  following  equation 

KeTe  =  Qe  (29) 

where 

Ke  =  Ke  -  gegl/Ke  (30) 

is  the  so-called  element  conductivity  matrix. 

Eq.  (29)  is  essentially  the  heat  conduction  equation  in  discrete  form  for  the  particle,  which  formally  defines 
the  discrete  thermal  element ,  and  all  the  associated  formulations  are  described  by  (19),  (20)  and  (28). 

A  very  important  feature  of  (29),  and  thus  this  new  thermal  element,  is  that  it  has  an  identical  form  as  a 
finite  element  thermal  element.  In  addition,  it  is  easy  to  verify  that 


Fig.  4.  Discrete  thermal  element. 
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Ke  =  KJe,  Kee  =  0  (31) 

i.e.  Ke  is  symmetric  and  singular,  as  in  the  case  of  finite  element  formulations.  Clearly  only  temperatures  are 
the  primary  variables,  while  the  local  fluxes  will  be  eliminated  in  the  subsequent  global  assembly  stage.  These 
features  are  different  from  the  standard  boundary  elements  where  both  nodal  temperatures  and  fluxes  are  the 
system  variables  and  the  resulting  matrices  are  unsymmetric  in  nature. 

Fig.  4  depicts  a  general  discrete  thermal  element,  where  the  number  of  nodes  n  equals  the  number  of  neigh¬ 
boring  particles  in  contact  and  thus  can  vary  from  element  to  element;  each  contact  zone  i  is  represented  by  a 
node  and  the  list  of  all  the  contact  zones  in  any  order  serves  as  the  element  connectivity;  while  the  contact  posi¬ 
tions  9i  and  angles  a*  are  the  geometric  properties  of  the  element. 

In  other  words,  each  circular  particle  can  be  viewed  as  a  special  shaped  finite  thermal  element,  although 
both  approaches  stem  from  different  principles.  The  immediate  consequence  is  that  the  subsequent  procedure 
in  this  discrete  thermal  element  method  to  model  the  assembled  particle  system  can  follow  the  same  steps  as 
used  in  FEM. 


3.3.  Thermal  resistance  matrix  He 


The  main  computational  issue  associated  with  the  formulation  of  the  discrete  thermal  element  proposed 
above  is  the  evaluation  of  the  coefficients  of  the  resistance  matrix,  htj.  Two  distinct  cases  need  to  be  considered: 


Case  1:  i  =  j.  In  this  case,  the  formula  for  ha  is 

e-(j) 


hii 


1  n  r 

4nKixf  F  J_a 


In 


sin- 


d(pdd 


(32) 


Clearly,  ha  is  solely  determined  by  the  contact  angle  a*.  As  the  integrand  In  |  sin (0  —  0)/2|  is  singular  at 
<fi  =  6,  some  caution  must  be  taken  when  numerically  evaluating  the  above  double  integral.  It  can 
however  be  simplified  into  a  single  integral  as  follows: 


2 

hu  = - /  (1  —  x)  ln[sin(az-x)]dx 

71K  J  o 


(33) 


Then  the  singularity  appears  only  at  x  =  0.  This  integral  may  be  effectively  solved  by  employing,  for 
instance,  a  special  Gaussian  rule  for  singular  functions  (see  [21]  for  more  detail).  Further  effort  has  led 
to  an  explicit  approximation  formula  for  hti  with  very  high  accuracy  as  follows: 


hu 

71K 


3  a/  d\  a \  a] 

'  11  ai  +  2  +  36  +  2700  +  79, 380  +  252, 000 


(34) 


Fig.  5a  shows  the  value  of  hu  as  a  function  of  the  contact  angle  a*  e  (0,  n/2]  (assuming  k  =  1),  where 
hu  is  monotonically  reduced  from  infinite  at  az-  =  0  to  the  minimum  value  of  0.356322  at  a*  =  n/2.  The 
higher  order  terms  of  a ,•  in  (34)  may  not  be  necessary  for  small  contact  angles.  Fig.  5b  depicts  the  accu¬ 
racy  of  the  above  formula  (34)  for  different  contact  angles,  along  with  the  accuracy  of  using  only  up  to 
the  terms  d\  and  a?.  Clearly  using  only  the  first  three  terms  (up  to  a?)  are  normally  sufficient  for  con¬ 
tact  angles  smaller  than  30°. 

Case  2:  i  ^  j.  In  this  case,  the  integrand  In  |  sin(A0Iy  +  6  —  </>)/2|  is  well-behaved  and  the  normal  2x2  Gauss¬ 
ian  quadrature  over  a  rectangular  domain  [—a*,  aj  x  [—ay,  ay]  can  be  used  to  compute  h tj  to  a  high  level 
of  accuracy.  For  small  contact  angles  a*  and  ay,  using  only  one  Gaussian  point  at  (0,0),  the  centre  of 
the  integration  domain,  may  be  sufficient.  Thus,  /zzy  can  be  evaluated  by  the  following  first  order 
approximation: 

hy  = — —  In  |  sin(0f  —  0j)/2\  (35) 


i.e.  hij  is  mainly  determined  by  the  angle  difference  between  the  two  contact  zones.  Given  that  the  two 
contact  zones  involved  have  the  same  contact  angle,  the  dependence  of  htj  (relative  to  hu)  on  the  angle 
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Fig.  5.  (a)  Value  of  hu  vs.  contact  angle;  (b)  accuracy  of  the  approximation  formula  (34). 


difference  \6t  —  6j\  is  illustrated  in  Fig.  6  for  three  different  contact  angles  a  =  2.5°,  5°  and  10°.  Smaller 
angle  difference,  i.e.  closer  contact  zones,  corresponds  to  a  larger  h tJ,  implying  a  stronger  coupling  be¬ 
tween  the  two  zones;  whereas  a  larger  angle  difference  results  in  a  weaker  coupling  effect.  When  the 
two  zones  are  separated  by  180°,  is  almost  zero,  is  generally  about  one  order  smaller  than  hu 
for  a  reasonably  small  contact  angle  as  the  angle  difference  between  any  two  contact  zones  is  normally 
larger  than  60°. 

Remark:  Formulae  (34)  and  (35)  not  only  demonstrate  that  the  thermal  resistance  matrix  of  the  element  can 
be  highly  effectively  evaluated,  but  also  identify  some  intrinsic  characteristics  of  heat  conduction  in  a  circular 
particle  -  the  thermal  resistance  properties  of  a  particle  can  be  solely  determined  by  its  thermal  contact  geo¬ 
metric  configuration,  in  which  the  contact  angles  a*  play  a  major  role,  and  the  contact  positions  6t  play  a  minor 
role,  while  the  size  of  the  particle  has  no  direct  effect. 
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Fig.  6.  Values  of  h^/hu  vs.  angle  differences. 


These  features  may  not  be  readily  apparent,  as  existing  models  often  take  particle  size  and  distance  between 
two  contacting  particles  as  the  main  parameters  to  determine  the  thermal  resistance/conductivity  of  a  particle. 
It  is  true  that  for  two  discrete  particles  brought  together  under  mechanical  forces,  the  corresponding  config¬ 
uration  of  the  contact  zone  is  closely  related  to  the  geometrical  properties  of  the  particles,  i.e.  the  contact  angle 
a  is  a  function  of  the  particle  sizes  and  their  distance.  However,  for  particles  that  are  ‘cemented’  or  bonded 
together,  the  distance  between  the  particles  is  irrelevant  to  the  contact  configuration.  Therefore,  the  use  of 
both  contact  angles  and  positions  in  this  approach  is  a  simple  and  universal  representation  of  the  thermal  con¬ 
tact  geometrical  features  of  a  particle. 

3.4.  The  solution  procedure  of  the  discrete  thermal  element  method 

The  discrete  thermal  element  method  shares  the  same  form  as  the  finite  element  formulation,  and  can  follow 
the  same  computational  steps  used  in  the  conventional  finite  element  procedure  to  solve  the  heat  conduction 
problem  in  a  particle  system.  The  main  steps  include:  (1)  preprocessing  to  obtain  all  element  configurations 
and  the  relevant  global  information;  (2)  evaluating  element  conductivity  matrices;  (3)  assembling  the  element 
contributions  to  form  a  global  linear  system  of  equations;  (4)  imposing  boundary  conditions;  (5)  solving  the 
equations  to  find  the  discrete  temperature  values;  and  (6)  postprocessing,  if  necessary,  to  compute  the  average 
temperature  for  each  particle  as  well  as  the  temperature  at  any  point. 

Based  on  the  positions  and  radii  of  particles,  the  contact  links  (zones)  between  the  particles  can  be  estab¬ 
lished  by  employing  any  search  algorithm  commonly  used  in  DEM  (see,  for  instance  [22]).  Each  link  is  equiv¬ 
alent  to  a  node  in  a  finite  element  model  and  is  thus  assigned  a  global  nodal  number.  The  total  number  of 
unknowns  in  the  global  system  is  equal  to  the  number  of  contact  links,  which  is  about  two  times  the  number 
of  particles.  The  list  of  the  contact  links  associated  with  a  particle  serves  as  its  element  connectivity,  and  the 
contact  positions  and  angles  are  the  associated  element  geometric  properties. 

Assuming  an  average  temperature  at  the  contact  zone  and  a  zero  net  flux,  i.e.  no  external  heat  source 
applied,  the  following  global  heat  conduction  equations  can  be  established  by  assembling  the  contributions 
from  all  the  particles: 

KT  =  Q  (36) 

where  K  is  the  global  conductivity  matrix;  T  is  the  global  temperature  vector  including  all  the  average  temper¬ 
ature  values  at  the  contact  regions;  and  Q  is  the  global  flux  (load)  vector  with  the  majority  of  the  components 
being  zero  if  internal  heat  sources  are  not  considered.  The  total  number  of  global  unknowns  equals  the  total 
number  of  contact  pairs  in  the  particle  system.  Note  that  K  is  symmetric  and  sparse. 
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With  appropriate  boundary  conditions  imposed,  the  resulting  global  system  of  equations  (36)  can  be  solved 
by  employing  one  of  the  advanced  sparse  matrix  solution  algorithms.  In  particular,  when  the  current  method 
is  coupled  with  DEM,  an  explicit  (pseudo-)time  integration  scheme,  similar  to  the  central  difference  algorithm, 
can  also  be  adopted.  It  has  been  established  that  under  certain  conditions  the  central  difference  time  integra¬ 
tion  is  equivalent  to  the  conjugate  gradient  method  [23].  Clearly,  systems  with  many  thousands  of  particles  can 
be  effectively  modelled. 

After  the  temperatures  at  the  contact  zones  are  obtained,  the  average  temperature  of  each  particle,  T0,  can 
be  readily  calculated  according  to  formula  (28).  The  element  flux  vector  Qe  can  be  computed  from  Eq.  (24). 
The  temperature  at  any  position  within  the  particle  domain  can  also  be  obtained  as  follows: 


T(r,d)  =  - 


1 

2kk 


In 


-  cos(0  -  6,  -<t>)  +  (£) 


d  <fi  +  T0 


(37) 


The  whole  solution  procedure  of  the  proposed  discrete  thermal  element  method  for  modelling  heat  condition 
in  a  particle  system  is  summarized  in  Box  1 . 


Box  1:  Discrete  thermal  element  method 

1.  Pre-processing 

•  Given  a  particle  assembly,  establish  the  contact  links  between  the  particles. 

•  For  each  particle,  identify  all  the  associated  contact  links/zones  and  determine  the  corresponding 
contact  positions  (0Z)  and  the  contact  angles  (a*). 

2.  Element  thermal  conductivity  computations.  Loop  over  all  particles, 

•  Evaluate  ha  and  h y  from  (34)  and  (20)  to  form  He. 

•  Compute  Ke  =  H”1. 

•  Determine  the  element  conductivity  matrix  Ke  based  on  (30). 

3.  Assemble.  Based  on  the  global  contact  links,  assemble  the  element  matrices  Ke  to  form  the  matrix  K. 

4.  Impose  boundary  conditions  to  (36). 

5.  Solve  the  global  system  of  equations  (36)  to  obtain  the  temperatures  at  the  contact  zones. 

6.  Post-processing,  if  necessary. 


4.  Additional  issues 

4.1.  Approximation  errors 

Unlike  finite  or  boundary  elements,  there  is  no  domain  discretization  involved  in  the  current  discrete  ther¬ 
mal  elements,  thus  no  discretization  error  will  occur.  Modelling  errors  however  stems  from  two  sources. 

Apart  from  the  slightly  inaccurate  geometry  description  for  the  contact  zones,  the  most  significant  assump¬ 
tion  made  in  the  development  of  the  discrete  thermal  element  formulations  is  that  the  flux  distribution  is  con¬ 
stant  over  the  contact  arc,  which  is  not  the  case  in  reality.  In  order  to  fully  understand  the  thermal  behavior 
within  the  contact  zone,  a  theoretical  analysis  is  conducted. 

Consider  two  contacting  disks  with  radii  RA  and  RB  and  the  corresponding  contact  angles  aA  and  aB,  respec¬ 
tively,  as  shown  in  Fig.  7a.  A  pair  of  point  heat  sources,  Q  and  —Q  are  applied  to  Points  E  and  F  at  the  far 
ends  of  the  two  disks,  respectively.  The  flux  distribution  along  the  arc  CHD  is  found  to  be 

gW  =  ~2«SaVV^')’  ^e(-aA’aA)  (38) 

where 

=  2«A^sin(«A)  sin(7t/?) _ +jf _ 

np  [^+^-2f^fYcos(2np)]{flf2)l-p 

with  p  =  n/y,  y  =  2n  -  aA  -  aB,  f  x  =  sin[(aA  +  </>)/ 2],  f2=  sin[(aA  -  4>)/2\. 


(39) 
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Fig.  7.  (a)  A  two-disk  contact  system;  (b)  comparison  of  flux  distribution  on  arc  CHD. 


This  solution  reveals  that  there  are  two  singular  points  at  </>  =  d=aA  where  the  flux  is  infinite. 

For  a  special  case  where  Ra  =  Rb  =  R,  aA  =  aB  =  a  =  7i/18  (=  10°)  and  Q  is  chosen  to  be  2otAR  so  that 
q((/))  =  —w((/)),  the  distribution  of  the  flux  along  the  arc  CHD,  computed  based  on  formula  (39),  is  depicted 
in  Fig.  7b.  Due  to  symmetry,  only  one  half  is  displayed.  To  assess  the  correctness  of  the  analytical  solution,  the 
FEM  result  obtained  by  using  a  very  fine  mesh  is  also  shown  in  the  figure.  An  excellent  agreement  between  the 
analytical  and  FEM  solutions  is  clearly  observed.  A  constant  flux  of  1,  corresponding  to  the  assumption  made 
in  the  earlier  development,  is  also  drawn  for  comparison. 

Fig.  7b  shows  that  the  actual  flux  is  fairly  constant  in  the  central  region  but  gradually  increases  towards  the 
end  of  the  arc  (the  intersection  point  of  the  two  disks),  at  which  the  flux  becomes  infinite.  A  small  discrepancy 
at  the  central  region  and  a  large  difference  towards  the  end  between  the  constant  flux  assumption  and  the 
actual  distribution  are  evident.  In  other  words,  formula  (32)  overestimates  ha  and  thus  introduces  solution 
errors.  Note  however  that  any  assumption  made  on  the  flux  distribution  will  have  little  effect  on  the  off-diag¬ 
onal  coefficients  p-,  as  long  as  the  two  contact  zones  concerned  are  not  too  close. 

Although  singularity  is  present  in  the  actual  flux  distribution,  the  temperature  distribution  in  the  contact 
zone  is  very  close  to  a  constant.  In  fact,  the  analytical  solution  indicates  that  for  two  disks  of  the  same  size, 
the  temperature  along  the  chord  CD  (refer  to  Fig.  7a)  must  be  constant. 

Based  on  the  above  analysis  and  observations,  two  solutions  can  improve  the  modelling  accuracy  of  the 
discrete  thermal  element  method.  In  the  first  approach,  the  analytical  flux  distribution  function  (38)  is  directly 
used  in  (17),  which  leads  to  the  following  new  formula  for  ha 


1  p  p 

hu  = - 2  /  /  w(^) ln 

mcy.f  J_a.  J_Xj 


.  6-4> 

sin — - — 


d(f>d0 


(40) 


where  w ((/)),  defined  in  (39),  can  be  viewed  as  a  weighting  function.  This  new  formula  should  provide  an  im¬ 
proved  value  for  p  in  principle.  However,  there  are  some  computational  implications  for  its  use  in  practice. 
Firstly,  it  is  complex  to  evaluate  w(</>).  This  issue  is  partially  resolved  for  small  contact  angles  by  using  a  sim¬ 
pler  approximate  function  to  replace  w(</>): 


_ 8/tei _ 

7i[(aA  +  4>)  ^  +  (aA  —  4>)  P]  (oc2a  —  </>2) 


(41) 


The  excellent  accuracy  of  this  approximation  is  demonstrated  in  Fig.  7b  for  the  special  conditions  given  ear¬ 
lier.  Secondly,  due  to  the  presence  of  the  additional  singularity  in  the  integral  (40),  extra  effort  has  to  be  made 
to  numerically  evaluate  the  integral.  Although  robust  numerical  procedures  can  be  derived  to  deal  with  the 
singularities  involved  in  the  integration,  a  highly  effective  formula  similar  to  (34)  may  no  longer  be  available. 
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For  the  above  reasons,  the  second  approach  for  improving  the  accuracy  of  the  discrete  thermal  element 
modelling  is  adopted.  It  is  proposed  that  instead  of  introducing  the  weighting  function  in  the  computation 
of  hti ,  a  correction  factor  X  is  directly  applied  to  reduce  the  value  of  ha  instead: 

ha  =  i)hu  (42) 

The  correction  factor  X  is  determined  by  directly  comparing  the  values  evaluated  by  the  two  formulae  (34)  and 
(40)  for  the  special  conditions  given  in  the  above,  with  the  latter  being  considered  as  the  exact  value.  Table  1 
lists  the  values  of  the  correction  factor  for  a  set  of  contact  angles.  The  corresponding  value  for  an  arbitrary 
contact  angle  can  be  obtained  by  a  simple  interpolation  of  the  values  given  in  the  Table.  Note  that  generally 
2-3%  of  error  is  present  under  the  special  conditions  in  the  original  DTEM. 

This  approach  retains  the  simplicity  and  effectiveness  of  the  original  method,  and  the  resulting  formula  is 
referred  to  as  the  enhanced  DTEM.  Note  that  the  correction  values  are  established  under  the  special  heating 
conditions,  therefore  the  approximation  error  of  the  DTEM  cannot  be  completely  eliminated  in  general  situ¬ 
ations.  The  performance  of  this  enhanced  version  will  be  assessed  through  numerical  validations  to  be  pre¬ 
sented  in  the  next  section. 

4.2.  Thermal  contact  resistance 

In  the  above  development,  a  perfect  thermal  contact  condition  is  assumed  between  two  contact  surfaces. 
However,  due  to  the  surface  roughness  as  well  as  other  complex  factors,  a  certain  degree  of  thermal  resistance 
is  always  present.  Such  a  phenomenon  can  be  readily  accommodated  in  the  current  discrete  thermal  element 
framework.  Basically  thermal  contact  resistance  can  be  simply  modelled  as  a  thermal  ‘bar’  or  'spring’.  When  a 
thermal  resistance  of  RTC  is  taken  into  account  at  the  contact  zones  of  two  particles,  an  equivalent  thermal  bar 
element  is  introduced  as 


'qF 
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in  which  and  T®  are,  respectively,  the  temperatures  at  the  contact  zones  of  the  two  particles;  and 
are  the  corresponding  heat  fluxes. 

This  thermal  contact  element  can  be  assembled  as  an  individual  element  into  the  global  conductivity  matrix 
in  the  normal  way.  The  only  difference  is  that  t\V)  and  T are  now  two  different  variables,  and  thus  the  total 
order  of  the  resulting  global  thermal  conduction  equations  will  be  doubled  if  thermal  resistance  is  considered 
at  every  contact  zone. 

5.  Numerical  validation  and  illustration 

Two  examples  are  provided  below  to  further  assess  the  solution  error  of  the  proposed  method  and  to  illus¬ 
trate  its  effectiveness. 

5.1.  Finite  element  validation 

The  solution  accuracy  of  the  proposed  discrete  thermal  element  method  with  the  two  versions  is  first 
assessed  against  the  finite  element  method.  An  assembly  of  12  particles  confined  within  an  unit  square  with 
two  vertical  walls  having  prescribed  temperatures,  as  shown  in  Fig.  8a,  is  used  as  an  example.  The  conductivity 
k  is  set  to  be  1.  The  dimensionless  coordinates  of  the  centres  and  radii  of  all  the  particles  are  listed  in  Table  2. 


Table  1 

Values  of  the  correction  factor  7  for  different  contact  angles  a 

a  (°)  1  2.5  5  7.5  10  12.5  15  17.5  20 

2  0.9799  0.9762  0.9729  0.9707  0.9693  0.9682  0.9675  0.9672  0.9671 
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For  geometrical  overlap,  there  are  in  total  17  contact  links/zones  between  the  particles,  plus  4  additional  con¬ 
tact  links  with  the  walls.  The  numberings  of  the  particles  and  the  contact  links  are  also  given  in  the  figure.  The 
thermal  conduction  between  the  particles  and  the  walls  can  be  treated  similarly  to  a  particle/particle  case.  The 
left  and  right  walls  have  a  prescribed  temperature  of  T\  =  1°  and  T2  =  0°,  respectively.  All  the  contact  angles, 
a,  are  within  a  narrow  range  of  5. 2-6. 2°.  The  values  of  n  vary  between  2  and  4;  and  the  total  number  of  the 
global  system  of  equations  is  17  after  imposing  the  boundary  conditions. 

The  corresponding  finite  element  model  using  triangular  elements  is  set  up  with  all  the  contact  zones  being 
represented  accurately.  Several  different  mesh  sizes,  ranging  from  0.02  to  0.0025,  are  used  to  provide  a  wider 
assessment  to  the  proposed  DTEM.  Fig.  8b  shows  the  coarsest  mesh  (mesh  size  =  0.02)  with  3839  nodes  and 
2141  elements.  For  this  mesh,  each  contact  zone  is  discretized  by  one  element.  The  results  obtained  from  the 
finest  mesh  (mesh  size  =  0.0025),  consisting  of  122,579  nodes  and  241,538  elements,  are  adopted  as  the  exact 
values.  Note  that  no  effort  has  been  made  to  use  mesh  adaptivity  to  optimize  the  finite  element  meshes.  Thus, 
the  sizes  of  the  FE  meshes  used  in  terms  of  elements  and  nodes  are  indicative  only. 

The  first  comparison  made  between  the  FE  models  and  the  two  discrete  thermal  element  versions,  standard 
and  enhanced,  is  for  the  computed  effective  thermal  conductivity,  k ,  of  the  assembly,  which  is  defined  as 

*  =  Q/{TX  —  To)  =  Q  (44) 

where  Q  is  the  total  heat  flux  flowing  into  the  assembly. 

For  the  enhanced  DTEM  version,  a  correction  factor  is  applied  to  each  diagonal  term  ha.  The  factor  for  the 
corresponding  contact  angle  is  obtained  by  linearly  interpolating  the  values  presented  in  Table  1. 

Table  3  summarizes  the  computed  effective  conductivity  for  the  different  models  together  with  the  errors 
relative  to  the  FE  model  of  mesh  size  =  0.0025.  It  can  be  seen  that  the  standard  DTEM  achieves  a  relative 
error  of  5.16%,  which  is  slightly  larger  than  that  of  the  finite  element  model  with  mesh  density  of  0.01  but 
smaller  than  that  of  0.02.  This  is  already  an  acceptable  level  for  most  engineering  problems.  It  is  also  shown 
that  the  enhanced  DTEM  further  reduces  the  error  by  more  than  half  to  2.37%  that  is  close  to  the  accuracy 
achieved  by  the  FE  model  with  mesh  size  =  0.005.  These  results  are  thus  very  encouraging,  particularly  when 
taking  into  consideration  the  simplicity  and  the  very  small  scale  of  the  DTEM  model.  The  remaining  error  is 
mainly  due  to  the  inaccurate  description  of  the  contact  geometry  which  will  be  investigated  later. 

The  average  temperatures  at  the  17  contact  zones  are  then  compared.  Table  4  lists  the  results  from  both 
versions  of  DTEM  and  FEM,  where  the  FEM  values  are  taken  from  the  finest  mesh,  and  DTEM-S  and 
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(a)  Particle  and  contact  link/zone  numberings  (b)  Mesh  size=0.02;  elements=3839; 

nodes=2414 


Fig.  8.  Example  1-12  particle  system:  (a)  discrete  particle  model;  and  (b)  finite  element  mesh. 
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Table  2 


Example  1  -  coordinates  of  centres  and  radii  of  particles 


Number 

X-coordinate 

7-coordinate 

Radius 

1 

0.120000 

0.120000 

0.1206 

2 

0.379230 

0.140000 

0.1407 

3 

0.638459 

0.120000 

0.1206 

4 

0.860000 

0.256087 

0.1407 

5 

0.210758 

0.363645 

0.1407 

6 

0.120000 

0.607291 

0.1206 

7 

0.544396 

0.340799 

0.1206 

8 

0.741856 

0.509942 

0.1407 

9 

0.408143 

0.562238 

0.1407 

10 

0.860000 

0.763796 

0.1407 

11 

0.601482 

0.736077 

0.1206 

12 

0.280799 

0.811604 

0.1407 

Table  3 

Example  1  -  comparison  of  computed  effective  conductivity  k 

Model 

FEM 

DTEM 

Mesh  size/version 

0.02 

0.01 

0.005 

0.0025 

Standard 

Enhanced 

Effective  k 

0.3062 

0.2931 

0.2860 

0.2824 

0.2678 

0.2757 

Error  (%) 

8.41 

3.79 

1.27 

- 

5.16 

2.37 

Table  4 

Example  1  -  comparison  of  average  temperatures  at  contact  zones 

Zone  number 

1 

2 

3 

4 

5 

6 

7 

8 

9 

FEM 

0.7216 

0.7678 

0.4964 

0.6679 

0.5217 

0.2923 

0.4340 

0.2228 

0.7807 

DTEM-S 

0.7208 

0.7674 

0.4965 

0.6668 

0.5218 

0.2917 

0.4337 

0.2223 

0.7803 

Error  (%) 

0.104 

0.056 

0.030 

0.176 

0.011 

0.220 

0.067 

0.233 

0.044 

DTEM-E 

0.7211 

0.7674 

0.4964 

0.6672 

0.5217 

0.2918 

0.4338 

0.2226 

0.7803 

Error  (%) 

0.072 

0.057 

0.007 

0.116 

0.017 

0.158 

0.049 

0.073 

0.050 

Zone  number 

10 

11 

12 

13 

14 

15 

16 

17 

Average 

FEM 

0.6426 

0.7719 

0.3604 

0.5056 

0.2128 

0.4648 

0.6353 

0.2657 

- 

DTEM-S 

0.6420 

0.7719 

0.3598 

0.5054 

0.2124 

0.4643 

0.6350 

0.2649 

- 

Error  (%) 

0.088 

0.001 

0.152 

0.037 

0.173 

0.126 

0.050 

0.307 

0.110 

DTEM-E 

0.6424 

0.7717 

0.3600 

0.5054 

0.2127 

0.4646 

0.6352 

0.2652 

- 

Error  (%) 

0.038 

0.034 

0.119 

0.034 

0.056 

0.056 

0.021 

0.172 

0.066 

DETM-E  denote,  respectively,  the  standard  and  enhanced  versions  of  the  DTEM.  It  shows  that  both  DETM- 
S  and  DTEM-E  have  achieved  much  higher  solution  accuracies  (comparing  with  the  effective  thermal  conduc¬ 
tivity),  with  averages  of  0.11%  and  0.066%,  respectively.  Again  the  enhanced  DTEM  is  about  two  times  more 
accurate  than  the  standard  DTEM.  The  fact  that  much  higher  accuracy  is  obtained  for  the  temperature  dis¬ 
tribution  than  the  flux  distribution  further  indicates  that  the  current  enhanced  DTEM  still  underestimates  the 
particle  conductivity  by  a  relative  constant  factor  and  thus  there  is  room  to  reduce  the  modelling  error  of  the 
DTEM  further. 

For  illustrative  purposes,  Fig.  9  depicts  the  average  temperature  distribution  of  the  particles  simulated  by 
the  DTEM-E  and  the  temperature  contour  plot  by  the  FEM. 
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(a)  Averaged  particle  temperature 


(b)  Temperature  contour 


Fig.  9.  Example  1  -  Particle  temperature  distribution  comparison:  (a)  DTEM;  (b)  FEM. 
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(a)  Discrete  thermal  element  model 

Fig.  10.  Example  2 
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(b)  Average  particle  temperature  distribution 

1084  particle  system. 


5.2.  Illustrative  example 

The  effectiveness  of  the  DTEM  is  further  demonstrated  with  this  second  example.  The  general  problem  set¬ 
ting  is  similar  to  the  previous  one  except  that  the  assembly  consists  of  1084  particles  that  are  randomly  gen¬ 
erated  by  the  advancing-front  packing  algorithm  [24].  The  sizes  of  particles  are  evenly  distributed  in  the  range 
[0.01, 0.02].  The  average  coordination  number  of  the  particles  is  3.88  and  there  are  2102  contact  links  in  total  - 
which  is  the  total  number  of  unknown  variables  in  the  DTEM  system.  The  corresponding  discrete  thermal 
element  model  is  depicted  in  Fig.  10a,  while  the  average  particle  temperature  distribution  simulated  by  the 
DTEM  is  illustrated  in  Fig.  10b,  from  which  a  realistic  result  is  apparent.  Note  that  much  larger  systems,  with 
tens  of  thousands  of  particles,  can  also  be  easily  simulated. 

6.  Concluding  remarks 

The  present  work  has  developed  a  novel  discrete  thermal  element  approach  that  not  only  provides  a  simple 
and  accurate  heat  conduction  model  for  particles  but  also  can  be  employed  to  effectively  simulate  heat  con- 
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duction  in  systems  comprising  a  large  number  of  particles.  The  method  exhibits  a  number  of  attractive 
features. 

Each  particle  is  treated  as  an  individual  element  with  the  number  of  (temperature)  unknowns  equal  to  the 
number  of  particles  that  it  is  in  contact  with.  The  element  thermal  conductivity  matrix  can  be  very  effectively 
evaluated  and  is  entirely  dependent  on  the  positions  and  contact  angles  of  the  associated  contact  zones/links, 
which  are  identified  as  the  main  characteristics  of  heat  conduction  in  circular  particles.  The  resulting  global 
system  of  equations  is  determined  by  the  topology  of  the  particle  contact  links  and  its  order  equals  the  number 
of  links. 

Both  the  element  and  global  conductivity  matrices  share  the  same  form  and  properties  (symmetry  and  spar¬ 
sity,  for  instance)  as  its  conventional  thermal  finite  element  counterpart.  Each  particle  can  thus  be  considered 
as  a  special  finite  element  and,  particularly,  the  whole  solution  procedure  follows  exactly  the  same  steps  as  that 
involved  in  the  finite  element  analysis.  Another  advantage  is  that  the  discrete  thermal  elements  can  be  com¬ 
bined  with  finite  elements  and  be  readily  incorporated  into  existing  finite  element  programs. 

There  is  no  discretization  error  associated  with  the  current  approach.  Apart  from  a  slightly  inaccurate 
geometry  description  for  the  contact  zones  for  small  contact  angles,  the  most  significant  assumption  made 
in  the  development  of  the  discrete  thermal  element  formulation  is  that  the  flux  is  constant  over  the  contact 
zone.  Based  on  the  theoretical  analysis,  a  simple  enhanced  version  has  been  proposed  to  deal  with  the  issue. 
The  numerical  assessment  conducted  indicates  that  the  standard  version  performs  exceptionally  well  for  the 
temperature  distribution,  exhibiting  about  5%  of  error  for  the  flux  for  the  given  example,  which  is  already 
within  the  acceptable  level  for  engineering  practice.  The  enhanced  version  further  improves  the  accuracy. 

Finally,  the  most  important  contribution  of  the  present  work  is  the  development  of  a  new  numerical  frame¬ 
work,  within  which  more  interesting  and  challenging  thermal  problems  in  particle  systems  can  be  more  effec¬ 
tively  tackled.  For  instance,  3D  spherical  particles  can  be  treated  in  a  similar  manner;  extension  to  transient 
problems  is  possible;  a  fully  coupled  mechanical-thermal  analysis  can  be  carried  out;  and  further  coupling 
with  fluid  flows  becomes  more  tractable.  These  problems  are  currently  being  pursued  and  will  be  reported 
later. 
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